Buckling analysis of nanobeams with exponentially varying stiffness by differential quadrature method
Chakraverty S, Behera Laxmi
Department of Mathematics, National Institute of Technology Rourkela, Odisha, India

 

† Corresponding author. E-mail: sne_chak@yahoo.com

Abstract

We present the application of differential quadrature (DQ) method for the buckling analysis of nanobeams with exponentially varying stiffness based on four different beam theories of Euler–Bernoulli, Timoshenko, Reddy, and Levison. The formulation is based on the nonlocal elasticity theory of Eringen. New results are presented for the guided and simply supported guided boundary conditions. Numerical results are obtained to investigate the effects of the nonlocal parameter, length-to-height ratio, boundary condition, and nonuniform parameter on the critical buckling load parameter. It is observed that the critical buckling load decreases with increase in the nonlocal parameter while the critical buckling load parameter increases with increase in the length-to-height ratio.

1. Introduction

Proper knowledge of mechanical properties is quite important in the design of various kinds of materials.[13] Due to their excellent physical, mechanical, and electrical properties,[4] nanostructures have attracted much attention among the scientists/researchers to develop innovatory applications in the field of nanomechanics. Proper understanding of their mechanical behavior is a key factor in the production of such engineering structures. Among these nanostructures, nanobeams attract more attention due to their great potential in engineering applications such as nanowires, nanoprobes, atomic force microscope (AFM), nanotube resonators,[5] nanoactuators,[6] and nanosensors. When nanostructure elements are subjected to compressive in-plane loads, these structures may buckle. Previously researchers have proposed that the surface stress itself can induce the buckling of a nanobeam with large deformation.[7] Figure 1 shows a schematic of buckling of a nanobeam with large deformation.

Fig. 1. (color online) Schematic of the buckling process of a nanobeam accompanied with surface stress variation.[7]

For a proper design of nanobeams, it is important to accurately predict the buckling characteristics of the nanobeams. The most challenging area in the field of nanomechanics is the study of structures at very small length scales. When the structures are of nanoscale size, classical continuum theories cannot accurately predict the mechanical behavior of such structures. Hence various size-dependent continuum theories such as strain gradient, couple stress, micropolar, and nonlocal elasticity have been developed to capture the size effects. Among these theories, the nonlocal elasticity theory proposed by Eringen[8] has gained much attention among the researchers because of its efficiency and simplicity to analyze the behavior of nanostructures. Applications of nonlocal continuum mechanics can be found in the areas of lattice dispersion of elastic waves, fracture mechanics, dislocation mechanics, wave propagation in composites, etc. It has been observed that the results obtained by the nonlocal continuum models are different from those previously obtained by the classical continuum theories, which indicates the influence of the size effects on the behavior of structures at nanoscale. Therefore, inclusion of the size effects in the buckling analysis of nanobeams is necessary. The presence of a nonlocal parameter in the constitutive relation indicates the inclusion of small scale effect. The nonlocal parameter has been calibrated using molecular dynamics to find its actual values.[9,10] The effect of the nonlocal parameter on the frequency and mode shapes is significant in the case of nanobeams.[11] Various relevant studies in the context of buckling analysis of nanobeams based on the nonlocal elasticity theory are shown in the subsequent paragraphs.

Buckling analysis of nano rods or tubes has been investigated analytically by Wang et al.[12] based on nonlocal Euler–Bernoulli and Timoshenko beam theories. Analytical results have been shown for four types of boundary conditions (simply-supported (SS), clamped-simply-supported (CS), clamped-clamped (CC), and cantilever (CF)). Reddy[13] reformulated various nonlocal beam theories such as Euler–Bernoulli beam theory (EBT), Timoshenko beam theory (TBT), Reddy beam theory (RBT), and Levison beam theory (LBT) based on the nonlocal elasticity theory and analytical results have been shown for SS nanobeams. Static nonlinear postbuckling response of nanobeams has been studied analytically by Emam[14] and the results show that as the nonlocal parameter increases, the critical buckling load reduces while the amplitude of buckling increases. Analytical methods have also been used in some of the studies[1520] related to nanobeams. One of the numerical methods, the Rayleigh–Ritz method, has been employed in most of the studies such as the buckling analysis of single-walled carbon nanotubes in thermal environments,[21] bending, buckling, and vibration problems of nonlocal Euler beams,[22] buckling analysis of single-walled carbon nanotubes subjected to different boundary conditions,[23] free vibration of nonhomogeneous Timoshenko nanobeams,[24] free vibration of rectangular nanoplates,[25] etc. A differential transformation method has been applied by Pradhan and Reddy[26] to predict the buckling behavior of a single-walled carbon nanotube (SWCNT) embedded on a winkler foundation. Sahmani and Ansari[27] developed a state-space method to study buckling of nanobeams based on Euler–Bernoulli, Timoshenko, and Levison beam theories. The effect of temperature on the buckling analysis of SWCNTs embedded in an elastic medium has been investigated by Narendar and Gopalakrishnan.[28] The nonlocal shell model has been used by Yan et al.[29] to study the nonlocal effect on axially compressed buckling of triple-walled carbon nanotubes under the influence of temperature. Some of the studies in context of buckling of single-layered graphene sheets based on plate theories were presented in Refs. [30]–[33].

Structural members with variable cross section are frequently used in civil, mechanical, and aeronautical engineering to satisfy architectural requirements. In practical cases such as space structures, this type of buckling analysis plays an important role in design. Many engineers currently design light slender members with variable cross sections to construct ever-stronger and ever-lighter structures. Unfortunately, design engineers are lacking proper knowledge on the design of nonuniform structural elements since most of the design specifications are available for uniform elements. Hence there is a need for buckling analysis of nonuniform structural elements. The structural design must ensure the stability of the initial equilibrium shape to evaluate correctly the critical loads corresponding to neutral and unstable equilibriums. Many structural elements have variable flexural rigidity which may be due to different reasons (technological, economical, etc.) and the present paper studies buckling analysis of such nonhomogeneous nanobeams.

Since conducting experiments with nanoscale size is quite difficult, so efforts are always being made by the researchers to develop efficient numerical or analytical methods for obtaining better results. The literature reveals that analytical methods have been used in most of the studies. However, it is not always possible to obtain analytical solutions for complicated geometries. Researchers have developed a few numerical methods in the buckling analysis of nanobeams. One must have sound knowledge about variational principles for applying numerical methods such as finite element method and Rayleigh–Ritz method. Again, subsequent application of variational principles often requires a proper understanding of principles of mechanics. This has motivated the search for an approximate computational technique for the direct solution of problems without recourse to variational principles. In this context, one may use the differential quadrature method.

The differential quadrature (DQ) method was introduced by Bellman and Casti[34] for the first time. This is an efficient numerical method for the solution of linear and nonlinear partial differential equations. Then Bert et al.[35] applied the DQ method in structural problems. Since then, researchers have been investigating linear and nonlinear structural problems using the DQ method. Different procedures have been used by the authors to implement the boundary conditions in the DQ method. Firstly, the δ technique was proposed by Bert et al.[36] to implement the boundary conditions. In this procedure,[37] one boundary condition is implemented at the boundary point and the other boundary condition is implemented at a distance δ from the boundary point. There are two major drawbacks to this approach. Firstly, since the implementation of the boundary condition at the δ point is an approximation of the true boundary condition which should be implemented at the boundary, accurate numerical solutions lie on the smaller value of δ. Secondly, smaller value of δ causes the solutions to oscillate since the weighting coefficient matrices become highly ill-conditioned. The above mentioned approach is suitable for clamped ends but not suitable for simply-supported and simply-supported-clamped. Seeing demerits of this approach, Bert[3842] proposed a new approach in applying the boundary conditions. In this approach, only one boundary condition is numerically implemented and the other boundary condition is built into the DQ weighting coefficient matrices. The following are some of the advantages of this approach. (i) The boundary conditions are properly satisfied since they are applied at the boundary points, therefore, better results may be obtained. (ii) The effect of δ on the results is eliminated. (iii) Excellent results are obtained with less computational effort.

On the other hand, the previous studies present investigations on the variable cross section in case of beams[43,44] and a few in case of nanostructures.[45] In case of buckling analysis of beams, the previous authors have used various ways of nonuniform cross section. In practical situations, stiffness of nanobeams may vary exponentially. To the best of the authors’ knowledge, till now, no work has been reported on buckling analysis of nanobeams having exponentially varying stiffness. It is also found that no literature on application of DQ in the buckling analysis of nanobeams subjected to different boundary conditions is available. In this article, the DQ method is employed to investigate nonuniform cross section on buckling analysis of nanobeams. Different beam theories such as EBT, TBT, RBT, and LBT are taken into consideration. Governing equations are converted into a single unknown function. Critical buckling loads are shown for various types of edge conditions such as SS,CS, CC, guided, and SS-guided boundary conditions. Effects of nonlocal parameter, nonuniformity, beam theories, boundary conditions, and length–height ratio have been investigated in this article.

2. Problem formulation

This study is carried out on the basis of four beam theories: Euler–Bernoulli beam theory, Timoshenko beam theory, Reddy beam theory, and Levinson beam theory. The present formulation is based on the nonlocal elasticity theory of Eringen.[8] Displacement fields, nonlocal constitutive equations, and stress–strain relations of the beam theories are well known and may be found in Ref. [13]. The governing equations which are converted to a single variable form are presented here. One should note that x, y, and z coordinates are taken along the length, width, and thickness of the beams, respectively (Fig. 2).

Fig. 2. (color online) Uniform beam with rectangular cross section and the coordinate system.
2.1. Euler–Bernoulli beam theory

The governing equation in terms of displacement may be written as[13]

where w is the transverse displacement, E is the Young’s modulus, I is the second moment of area, and is the applied axial compressive force. It may be noted here that is the nonlocal parameter, where and a denote the material constant and the internal characteristic length, respectively.

2.2. Timoshenko beam theory

The governing equations for TBT nanobeams are given by[13]

where φ is the rotation of the cross-section, G is the shear modulus, A is the cross sectional area, and is the shear correction factor. By eliminating φ from Eqs. (2) and (3), the governing equations can be transformed into an uncoupled differential equation as follows:

2.3. Reddy beam theory

In this case, the governing equations in terms of displacements may be written as[13]

where are the second, the fourth, and the sixth order moments of area about the y axis and are defined as , , , , , . Here and , where h is the height of the beam. one may note that higher order inertias, i.e., , and , are ignored in the formulation.[13]

By eliminating φ from Eqs. (5) and (6), the governing equations may be obtained in terms of displacement as

2.4. Levinson beam theory

The governing differential equations may be written as[13]

By eliminating φ from Eqs. (8) and (9), the uncoupled governing differential equation is
The notations in the above equations are already defined in the previous subsections.

Next, we assume an exponential variation of the flexural stiffness (EI) since the flexural stiffness of the SWCNT may not be constant for a geometrically non-uniform beam model of SWCNT. For this, we consider a beam with a nonuniform variation of the cross section along the length. Based on the proposed exponential variation, the flexural stiffness is defined as follows:

where is the mass moment of area at the left end and η is a positive constant.

3. Solution methodology

In this analysis, we have taken the following non-dimensional terms:

Below we show the non-dimentionalized forms of the governing differential equations for EBT, TBT, RBT, and LBT respectively:

Based on the above equations, the solution of the governing equations will be sought using the DQ method. A brief introduction of the method is given below.

3.1. Differential quadrature method

The derivatives of displacement function at a given discrete point i are approximated as

where i = 1, 2, …, N and N is the number of discrete grid points. Here , and are the weighting coefficients of the first, the second, the third, and the fourth derivatives, respectively. One should note that other higher order derivatives present in the governing equations are approximated in a similar way.

Determination of weighting coefficients Computation of weighting coefficient matrix is the key step in the DQ method. In the present investigation, we have used Quan and Chang’s[46] approach to compute weighting coefficients . Matrix may be computed by the following procedure: for ,

for ,
Once the weighting coefficients of the first order derivatives are computed, the weighting coefficients of the higher order derivatives can be obtained by simply matrix multiplication as follows:

Selection of mesh point distribution We assume that the domain is divided into intervals with coordinates of the grid points given as . These ’s have been computed by using Chebyshev–Gauss–Lobatto grid points, that is,

One may note that nonuniform grid points have been used.

Application of boundary condition The above matrices , and are converted into modified weighting coefficient matrices , and as per the boundary conditions which will be shown below. Five boundary conditions, i.e., SS, CS, CC, guided, and SS-guided boundary conditions, are considered in the present analysis.

Let us now denote

In view of the above, we now illustrate the procedure to find the modified weighting coefficient matrices below.

For the simply supported boundary condition, equation (15) may be rewritten in matrix form as

Firstly, by applying boundary condition , equation (21) becomes
or
For the second derivative, one has
By using the above relation for
where

Now, since , the third order derivative of the displacement function may be obtained as follows:

Here

Similarly, for the fourth order derivative,

where or

For the clamped-simply supported boundary condition, proceeding in a similar fashion as before, we have , , with , , with , , with .

For the clamped-clamped boundary condition, we have , , with , , with , , with .

For the guided boundary condition, we have , , with , , with , with

For the simply supported-guided boundary condition, we have , , with , , with , , with

It is noted that when substituting the derivatives in the governing differential equations, one has to use , , , as per the boundary condition. Substituting Eq. (15) into Eqs. (11)–(14), one may obtain the generalised eigen value problem as

where is the stiffness matrix and is the buckling matrix.

4. Numerical results and discussion

In this section, numerical results have been obtained by solving Eq. (23) using Matlab program developed by the authors. The DQ method has been employed and the boundary conditions are implemented in the coefficient matrix. It may be noted that the following parameters are taken for the computation: E = 70 Tpa, , h = 1, .

4.1. Validation

First, to validate the present method, we compare our results of the critical buckling load parameter with those available in the literature. For the validation, we consider uniform () nanobeams. To compare our results of EBT and TBT nanobeams with that of Wang et al.,[12] we consider a beam of diameter d = 1 nm, Young’s modulus E = 1 Tpa, and Poisson’s ratio . The comparison is shown in Table 1 for three types of boundary conditions (SS, CS, and CC). In this table, the critical buckling load parameter (in nN) for EBT and TBT nanobeams with is presented for various scale coefficients (0, 0.5 nm, 1 nm). Similarly, the critical buckling loads of RBT and LBT nanobeams are compared respectively with those from Refs. [14] and [27] in Table 2. For this purpose, the same parameters as those used in Refs. [14] and [27] are taken. It may be noted that the comparison for the RBT nanobeams has been made with length-to-height ratio and nonlocal parameter as 0, 1 nm, 2 nm, 3 nm while comparison for the LBT nanobeams has been done with and as 0, 0.5 nm, 1 nm 1.5 nm. In this table, we have considered SS and CC edge conditions. It is seen that the critical buckling load parameter decreases with the increase of the nonlocal parameter. From these tables, one may observe close agreement of results with those available in the literature.

Table 1.

Comparison of critical buckling load parameter (in units of nN) for EBT and TBT nanobeams.

.
Table 2.

Comparison of critical buckling load parameter for RBT and LBT nanobeams.

.
4.2. Convergence

To find the minimum number of grid points for obtaining accurate results, a convergence study is carried out for EBT and TBT nanobeams. To show how the solution is affected by the grid points, the variation of the critical buckling load parameter with the number of terms (N) is shown in Fig. 3. In this figure, we have taken nm, nonuniform parameter , and L = 10 nm. The convergence is shown for the simply supported edge condition only. From this figure, one may observe that convergence is achieved as we increase the number of terms. It may be noted that eleven grid points are sufficient to compute the results.

Fig. 3. (color online) Variation of with grid points.
4.3. Small scale effect

In this subsection, the significance of the scale coefficient is highlighted. First we define the buckling load ratio as calculated using nonlocal theory and / calculated using local theory. This buckling load ratio serves as an index to estimate quantitatively the small scale effect on the buckling solution. To state the importance of the scale coefficient, the variation of the buckling load ratio with scale coefficient is shown in Figs. 47. In these figures, we have considered RBT and LBT nanobeams with edge conditions as the guided and simply supported-guided boundary conditions and nonuniform parameter as . The results for different length-to-height ratios are shown. It is noticed from the figures that the buckling load ratios are less than unity. This implies that the application of the local beam model for the buckling analysis of carbon nanotubes would lead to overprediction of the buckling load if the small length scale effect between the individual carbon atoms is ignored. As the scale coefficient increases, the buckling loads obtained by the nonlocal beam model become smaller than those of its local counterpart. In other words, the buckling load parameter obtained by the local beam theory is larger than that obtained by the nonlocal beam theory. So the presence of the nonlocal parameter in the constitutive equation is significant in the field of nanomechanics. It is also observed that the small scale effect is affected by . This observation is explained as follows. When increases, the buckling load ratio comes closer to one. This implies that the buckling load parameter obtained by the nonlocal beam model comes closer to that furnished by the local beam model. Hence the small scale effect is negligible for a very slender carbon nanotube (CNT) while it is significant for short carbon nanotubes (CNTs). This implies that if we compare magnitude of small scale effect with length of the slender tube, the small scale (internal characteristic length) is so small that it can be regarded as zero.

Fig. 4. (color online) Variation of buckling load ratio with (RBT guided).
Fig. 5. (color online) Variation of buckling load ratio with (RBT SS guided).
Fig. 6. (color online) Variation of buckling load ratio with (LBT guided).
Fig. 7. (color online) Variation of buckling load ratio with (LBT SS guided).
4.4. Effect of nonuniform parameter

In this subsection, the effect of non-uniformity η on the critical buckling load parameter is illustrated. It may be noted that η takes positive values only and different value of η gives different cross section of the nanobeams. This analysis will help design engineers to have an idea of the values of the critical buckling load parameter. Figure 8 shows the variation of critical buckling load parameter with nonuniform parameter η. In this graph, we have considered EBT nanobeams with CS edge conditions having L = 50 nm. The results are shown for different scale coefficients. It is observed that with increase in the nonuniform parameter, the critical buckling load parameter decreases. This decrease is attributed to the decrease of the stiffness of the beam. It is also observed that with increase in the nonlocal parameter, the critical buckling load parameter decreases. This shows that the local beam theory () over predicts the buckling load parameter. Hence for better predictions of the buckling load parameter, one should consider the nonlocal theory. Next, to investigate the influence of the nonuniform parameter on the higher buckling mode, the variation of buckling load parameter with nonuniform parameter η is presented in Fig. 9. Here we have considered LBT nanobeams with SS edge condition having nm and L = 15 nm. It is seen from the figure that the buckling load parameter decreases with increase in the nonuniform parameter and this decrease is more significant in case of higher modes.

Fig. 8. (color online) Variation of with nonuniform parameter.
Fig. 9. (color online) Variation of with nonuniform parameter.
4.5. Effect of length-to-height ratio

Another important factor that design engineers should keep in mind is the length-to-height ratio. To investigate the effect of length-to-height ratio on the critical buckling load parameter, the variation of the critical buckling load parameter with is shown in Figs. 1013. The results are shown for different scale coefficients (0.5 nm, 1 nm, 1.5 nm). In these graphs, we have considered EBT and TBT nanobeams with edge conditions such as CS and CC. The numerical values have been obtained by taking . It is observed that the critical buckling load parameter increases with increase in . It is also noticed that the critical buckling load parameter decreases with increase in the scale coefficient. Hence one should incorporate the nonlocal theory in the buckling analysis of nanobeams. One may also notice that for a particular , the results obtained by TBT nanobeams are less compared to EBT nanobeams. This is due to the absence of transverse shear stress and strain in the EBT nanobeams. One may say that TBT nanobeams can be better predicted of the critical buckling load parameter than EBT nanobeams.

Fig. 10. (color online) Variation of with (EBT CC).
Fig. 11. (color online) Variation of with (EBT CS).
Fig. 12. (color online) Variation of with (TBT CC).
Fig. 13. (color online) Variation of with (TBT CS).
4.6. Effect of various beam theories

Modeling of nanostructures based on beam theories is one of the important areas in the field of nanotechnology. Various beam theories used in this investigation are discussed below. The first and simplest beam theory is EBT. In this theory, the transverse shear and transverse normal strains are ignored. The next theory developed in the hierarchy of beam theories is TBT. According to this theory, a constant state of transverse shear strain with respect to the thickness coordinate is assumed. Therefore, a shear correction factor is required in order to compensate the error due to the above assumption. In RBT and LBT, the transverse shear strain varies quadratically and also vanishes on the top and bottom planes of the beams. Thus the shear correction factor is not taken into consideration. EBT and TBT are the first order beam theories, while RBT and LBT are the third order beam theories. To investigate the effect of various beam theories such as EBT, TBT, RBT, and LBT on the buckling load parameter, the variation of the critical buckling load parameter with scale coefficient is shown in Fig. 14 for various types of beam theories. In this figure, the CS boundary condition is taken into consideration with L = 10 nm and . It is seen from the figure that EBT predicts higher critical buckling load parameter than the other types of beam theories. It is due to the fact that in EBT, the transverse shear and transverse normal strains are not considered. It is also noted that the beam theories such as TBT, RBT and LBT predict approximately the same results.

Fig. 14. (color online) Variation of with small scale coefficient.
4.7. Effect of boundary condition

For designing engineering structures, one must have proper knowledge about the boundary conditions. It will help engineers to have an idea without carrying out detailed investigation. Therefore, analysis of the boundary conditions is quite important. In this subsection, we consider the effect of the boundary condition on the critical buckling load parameter. Figure 15 depicts the variation of the critical buckling load parameter of the TBT nanobeam with the scale coefficient for different boundary conditions. This graph is plotted with L = 10 nm and . It is observed form the figure that CC nanobeams have the highest critical buckling load parameter and SS-guided nanobeams have the lowest critical buckling load parameter. It is also seen that SS and guided boundary conditions give almost the same results.

Fig. 15. (color online) Variation of with small scale coefficient.
5. Conclusion

Buckling analysis of nonuniform nanobeams based on four different beam theories including EBT, TBT, RBT, and LBT has been carried out. The nonuniform cross section is assumed by taking exponentially varying stiffness. The differential quadrature method has been employed and the boundary conditions have been substituted in the coefficient matrix. New results have been shown for two types of boundary conditions, i.e., guided and SS-guided boundary conditions. The numerical results are presented to show the effects of the nonuniform parameter, the nonlocal parameter, , the boundary condition, and the beam theories on the critical buckling load parameter. It is found that the critical buckling load parameter decreases with increase in the nonlocal parameter. It is also observed that the application of the Euler–Bernoulli beam theory would lead to overprediction of the buckling loads. One may also notice that CC nanobeams have higher buckling loads than the other types of boundary conditions. We conclude that for better prediction of buckling loads, one should incorporate the nonlocal theory.

Reference
[1]YangLHuLZhangX G2015Acta Phys. Sin.640134502
[2]ZhangWHuLZhangG2016Acta Phys. Sin.65024502
[3]DaiZ DLiS LLiD LZhuA J 2007 Chin. Phys. Lett. 24 1429
[4]RuudJ AJervisT RSpaepanF 1994 J. Appl. Phys. 75 4969
[5]PengHChangCAloniSYuzvinskyTZettlA 2006 Phys. Rev. Lett. 97 087203
[6]DubeyASharmaGMavroidisCTomassoneMNikitczukKYarmushM 2004 Journal of Computational and Theoretical Nanoscience 1 18
[7]LiuJSunJZuoP 2016 Acta Mechanica Solida Sinica 29 192
[8]EringenA C1972International Journal of Engineering Science101
[9]DuanWWangC MZhangY 2007 J. Appl. Phys. 101 024305
[10]HuangLHanQLiangY2012Nano7125003
[11]XuM 2006 Proc. Roy. Soc. A: Math. Phys. Eng. Sci. 462 2977
[12]WangC MZhangY YRameshS SKitipornchaiS 2006 J. Phys.D: Appl. Phys. 39 3904
[13]ReddyJ2007International Journal of Engineering Science45307
[14]EmamS A 2013 Applied Mathematical Modelling 37 6929
[15]AydogduM 2009 Physica 41 1651
[16]ThaiH T 2012 International Journal of Engineering Science 52 56
[17]ReddyJPangS 2008 J. Appl. Phys. 103 023511
[18]KumarDHeinrichCWaasA M 2008 J. Appl. Phys. 103 073521
[19]SimşekMYurtcuH 2013 Composite Structures 97 378
[20]YangYLimC W 2011 Nano 6 363
[21]AnsariRSahmaniSRouhiH 2011 Computational Materials Science 50 3050
[22]GhannadpourS A MMohammadiBFazilatiJ 2013 Composite Structures 96 584
[23]AnsariRSahmaniSRouhiH 2011 Phys. Lett. 375 1255
[24]BeheraLChakravertyS 2014 Meccanica 49 51
[25]ChakravertySBeheraL 2014 Physica 56 357
[26]PradhanSReddyG 2011 Computational Materials Science 50 1052
[27]SahmaniSAnsariR2011Journal of Mechanical Science and Technology252365
[28]NarendarSGopalakrishnanS 2011 Physica 43 1185
[29]YanYWangWZhangL 2010 Applied Mathematical Modelling 34 3422
[30]PradhanS C 2009 Phys. Lett. 373 4182
[31]PradhanS CMurmuT 2010 Physica 42 1293
[32]PradhanS C 2012 Sadhana 37 461
[33]NarendarS 2011 Composite Structures 93 3093
[34]BellmanRCastiJ 1971 Journal of Mathematical Analysis and Applications 34 235
[35]BertC WJangS KStrizA G 1988 AIAA J. 26 612
[36]JangS KBertC WStrizA G 1989 International Journal for Numerical Methods in Engineering 28 561
[37]ShuC2000Differential Guadrature and its Application in Engineering New YorkSpringer
[38]WangXBertC W 1993 Journal of Sound and Vibration 162 566
[39]BertCWangXStrizA 1994 Acta Mechanica 102 11
[40]WangXBertCStrizA 1993 Computers and Structures 48 473
[41]BertC WXinweiWStrizA G1993International Journal of Solids and Structures301737
[42]BertC WMalikM 1996 International Journal of Mechanical Sciences 38 589
[43]KonstantakopoulosT GRaftoyiannisI GMichaltsosG T2012Open Construction and Building Technology Journal61
[44]SapountzakisETsiatasG 2007 Engineering structures 29 675
[45]FarajpourADaneshMMohammadiM2007Physica E44719
[46]QuanJChangC 1989 Computers and Chemical Engineering 13 779